
library(ggplot2)
library(ggpubr)
library(haven)
library(dplyr)


folder = ""

inputs = paste0(folder, "\\2_Intermediate")

event_studies = read_dta(paste0(inputs, "\\Event Studies.dta"))
dataset = read_dta(paste0(inputs, "\\CC NY Comparison.dta"))
boxplots = read_dta(paste0(inputs, "\\Boxplots.dta"))

save_dir = paste0(folder, "\\3_Outputs\\")


################################################################################
# Plot Settings

labs_font = 25
axis_name_size = 18
sub_title_font = 22
legend_font = labs_font
tick_len = .35

custom_theme = theme(
  axis.title.x = element_text(size = labs_font),
  axis.text = element_text(size = axis_name_size),
  axis.title.y = element_text(size = labs_font),
  plot.title = element_text(size = sub_title_font),
  legend.text = element_text(size = legend_font),
  axis.ticks.length=unit(tick_len, "cm"),
  legend.direction = "horizontal",
  legend.position = "bottom",
  strip.text.x = element_text(size = labs_font),
  strip.background = element_blank(),
  legend.key.size = unit(4,"line"))


################################################################################
#
## COMBINED EVENT STUDIES, AND SINGLE-CITY EVENT STUDIES
## INCLUDES FIGURES 2, 3, 4, AND SUPPLEMENTARY FIGURES C.1, E.1
# 
################################################################################

# Graph-specific Settings
outcomes = c('Number of Hospital Visits', 'Number of Emergency Visits', 'Number of Mental Health Visits', 
             'Observed this Year', 'Any Payday Inquiry * 100', 'Number of Payday Inquiries','Any Payday Account * 100', 'Number of Payday Loans',
             'Earnings_comb', 'Employment_comb', 'Not at Eviction Address_comb', 'Any Homelessness Services_comb', 'Emergency Shelter Use_comb', 'Neighborhood Poverty_comb',
             'Credit Score_comb', 'No Open Revolving Account_comb', 'Any Auto Loan or Lease_comb', 
             'Total Balance: Col. and Del._comb','Financial Health Index_comb') 


outcome_labels = c('Number of Hospital Visits', 'Number of Emergency Visits', 'Number of Mental Health Visits', 
                   'Observed this Year', 'Any Payday Inquiry * 100', 'Number of Payday Inquiries','Any Payday Loan * 100', 'Number of Payday Loans',
                   'Earnings', 'Employment', 'Not at Eviction Address', 'Any Homelessness Services', 'Emergency Shelter Use', 'Neighborhood Poverty',
                   'Credit Score', 'No Open Revolving Account', 'Any Auto Loan or Lease', 
                   'Total Balance: Col. and Del.','Financial Health Index') 


save_names = c('Figure_4_Hospital', 'Figure_4_Emergency', 'Figure_4_Mental',
               'Supplementary data_Figure_E1_Obs this year', 
               'Supplementary data_Figure_C1_Payday Inq','Supplementary data_Figure_C1_Num Payday Inq', 
               'Supplementary data_Figure_C1_Payday Acc', 'Supplementary data_Figure_C1_Num Payday Loan',
               'Figure_2_Comb Earnings', 'Figure_2_Comb Employment', 'Figure_2_Comb Residential Mobility', 'Figure_2_Comb Homelessness', 
               'Figure_2_Comb Emergency', 'Figure_2_Comb Pov',
               'Figure_3_Comb Credit', 'Figure_3_Comb Zero Revolving', 'Figure_3_Comb Auto Loan', 
               'Figure_3_Comb Deli','Figure_3_Comb Index')


relative_quarter = c(1,1,1,
                     0,1,1,1,1,
                     1,1,0,0,0,0,
                     1,1,1,1,1)


all_graphs = list()

for (o in 1:length(outcomes)) {
  
  print(outcomes[o])
  
  df = event_studies %>% filter(outcome==outcomes[o])
  
  if (relative_quarter[o]) {
    x_breaks = seq(from = -16, to = 24, by = 4)
  }
  else {
    x_breaks = min(df$time):max(df$time)
  }
  
  x_labs = x_breaks
  
  if (outcomes[o]=='Number of Hospital Visits') {ylabels = c(0.16, 0.18, 0.20, 0.22, 0.24)}
  else if (outcomes[o]=='Number of Emergency Visits') {ylabels = c(0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20)}
  else if (outcomes[o]=='Number of Mental Health Visits') {ylabels = c(0.005, 0.010, 0.015, 0.020, 0.025)}
  
  else if (outcomes[o]=='Observed this Year') {ylabels = c(0.825, 0.850, 0.875, 0.900, 0.925, 0.950)}
  
  else if (outcomes[o]=='Any Payday Inquiry * 100') {ylabels = c(2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0)}
  else if (outcomes[o]=='Number of Payday Inquiries') {ylabels = c(0.075, 0.100, 0.125, 0.150, 0.175, 0.200)}
  else if (outcomes[o]=='Any Payday Account * 100') {ylabels = c(0.0, 0.1, 0.2, 0.3, 0.4, 0.5)}
  else if (outcomes[o]=='Number of Payday Loans') {ylabels = c(0.001, 0.002, 0.003, 0.004, 0.005, 0.006)}
  
  else if (outcomes[o]=='Earnings_comb') {ylabels = c(3200, 3400, 3600, 3800, 4000, 4200, 4400)}
  else if (outcomes[o]=='Employment_comb') {ylabels = c(0.50, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57)}
  
  else if (outcomes[o]=='Not at Eviction Address_comb') {ylabels = c(0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)}
  else if (outcomes[o]=='Emergency Shelter_comb') {ylabels = c(0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8)}
  else if (outcomes[o]=='Any Homelessness Services_comb') {ylabels = c(0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08)}
  else if (outcomes[o]=="Neighborhood Poverty_comb") {ylabels = c(0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27)}
  
  else if (outcomes[o]=='Credit Score_comb') {ylabels = c(530, 535, 540, 545, 550, 555, 560)}
  else if (outcomes[o]=='No Open Revolving Account_comb') {ylabels = c(0.42, 0.44, 0.46, 0.48, 0.50, 0.52, 0.54)}
  else if (outcomes[o]=='Total Balance: Col. and Del._comb') {ylabels = c(2000, 2200, 2400, 2600, 2800, 3000)}
  else if (outcomes[o]=='Any Auto Loan or Lease_comb') {ylabels = c(0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18)}
  else if (outcomes[o]=='Financial Health Index_comb') {ylabels = c(-0.30, -0.20, -0.10, 0.00, 0.10)}
  
  ylimits = c(min(ylabels), max(ylabels))
  
  print(outcomes[o])
  
  g = ggplot(df) + 
    theme_classic() +
    geom_vline(xintercept = 0, color = "black", alpha = .75, linetype = 'dashed') +
    geom_line(aes(x = time, y = evicted,  color = "Evicted", linetype = "Evicted"), size = 2)  +
    geom_line(aes(x = time, y = non_evicted, color = "Not Evicted", linetype = "Not Evicted"), size = 2) +
    scale_color_manual(name = "", 
                       values = c("Evicted" = "#1B9E77", "Not Evicted" = "#D95F02")) +
    scale_linetype_manual(name = "",
                          values = c("Evicted" = "solid", "Not Evicted" = "dashed")) +
    custom_theme + 
    scale_x_continuous(breaks = x_breaks, labels = x_labs) +
    scale_y_continuous(limits = ylimits, breaks = ylabels, labels = ylabels) +
    labs(x=ifelse(relative_quarter[o], 'Relative Quarter', 'Relative Year') , y=outcome_labels[o]) + 
    theme(legend.position="bottom")
  
  ggsave(paste0(save_dir, save_names[o], ".png"), plot=g, width=10, height=7)
  
  print(outcomes[o])
}



################################################################################
#
## MAKE FIGURE 6
#
################################################################################

lower_lim1 <- -0.2
upper_lim1 <- 0.2

lower_lim2 <- -0.2 
upper_lim2 <- 0.2


dataset_to_graph <- dataset
dataset_to_graph$Outcome <- factor(dataset_to_graph$Outcome,
                                   levels = c("Not at eviction address",
                                              "Emergency shelter",
                                              "Any homelessness services",
                                              "Neighborhood poverty rate",
                                              "Earnings",
                                              "Employment",
                                              "Financial health index",
                                              "Credit score",
                                              "No open revolving account",
                                              "Total balance: collections and delinquencies",
                                              "Any auto loan or lease"))

dataset_to_graph$Group <- factor(dataset_to_graph$Group,
                                 levels=c("Housing", "Labor", "Financial"))

iv_only <- dataset_to_graph %>% filter(Specification=="IV")

g <- ggplot(iv_only, aes(valuen, valuec, shape = Group, color = Outcome)) + 
  geom_point(size = 4, alpha = .9, color = "black") +
  geom_point(size = 3, alpha = .99) +
  geom_abline(intercept = 0, slope = 1, alpha = .5) +
  ylab("Cook County Estimates (standardized)") +
  xlab("New York Estimates (standardized)") +
  theme_bw() +
  xlim(c(lower_lim1, upper_lim1)) + ylim(c(lower_lim1, upper_lim1)) +
  scale_color_brewer(palette = "Spectral")

ggsave(paste0(save_dir, "Figure_6.png"), plot=g, height = 5, width = 8.7)




################################################################################
#
## COOK VS. NY EVENT STUDIES 
## INCLUDES SUPPLEMENTARY FIGURES E.2 AND E.3
# 
################################################################################

# Graph-specific Settings
outcomes = c('Earnings', 'Employment', 'Not at Eviction Address', 'Any Homelessness Services', 'Emergency Shelter Use', 'Neighborhood Poverty',
             'Credit Score', 'No Open Revolving Account', 'Any Auto Loan or Lease',
             'Total Balance: Col. and Del.', 'Financial Health Index' )

save_names = c('Supplementary data_Figure_E2_Earnings', 'Supplementary data_Figure_E2_Employment', 'Supplementary data_Figure_E2_Residential Mobility',
               'Supplementary data_Figure_E2_Homelessness', 'Supplementary data_Figure_E2_Emergency Shelter Use', 'Supplementary data_Figure_E2_Pov', 
               'Supplementary data_Figure_E3_Credit', 'Supplementary data_Figure_E3_Zero Revolving',
               'Supplementary data_Figure_E3_Auto Loan', 'Supplementary data_Figure_E3_Deli', 'Supplementary data_Figure_E3_Index')

relative_quarter = c(1, 1, 0, 0, 0, 0,
                     1, 1, 1, 
                     1, 1)


all_graphs = list()

for (o in 1:length(outcomes)) {
  
  print(outcomes[o])
  
  df = event_studies %>% filter(outcome==outcomes[o])
  df_cc = df %>% filter(geo=='cc')
  df_ny = df %>% filter(geo=='ny')
  
  if (relative_quarter[o]) {
    x_breaks = seq(from = -16, to = 24, by = 4)
  }
  else {
    x_breaks = min(df$time):max(df$time)
  }
  
  x_labs = x_breaks
  
  
  if (outcomes[o]=='Earnings') {
    cc_ylabels = c(3600, 3800, 4000, 4200, 4400, 4600, 4800, 5000, 5200)
    ny_ylabels = c(2800, 2900, 3000, 3100, 3200, 3300, 3400, 3500)
  }
  else if (outcomes[o]=='Employment') {
    cc_ylabels = c(0.56, 0.58, 0.60, 0.62, 0.64, 0.66)
    ny_ylabels = c(0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48)
  }
  else if (outcomes[o]=='Not at Eviction Address') {
    cc_ylabels = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)
    ny_ylabels = c(0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8)
  }
  else if (outcomes[o]=='Emergency Shelter Use') {
    cc_ylabels = c(0.00, 0.005, 0.01, 0.015, 0.02)
    ny_ylabels = c(0.00, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12)
  }
  else if (outcomes[o]=='Any Homelessness Services') {
    cc_ylabels = c(0.00, 0.01, 0.02, 0.03, 0.04)
    ny_ylabels = c(0.00, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12)
  }
  else if (outcomes[o]=='Neighborhood Poverty') {
    cc_ylabels = c(0.14, 0.16, 0.18, 0.2, 0.22)
    ny_ylabels = c(0.26, 0.27, 0.28, 0.29, 0.3)
  }
  else if (outcomes[o]=='Credit Score') {
    cc_ylabels = c(510, 520, 530, 540, 550)
    ny_ylabels = c(540, 550, 560, 570, 580)
  }
  else if (outcomes[o]=='No Open Revolving Account') {
    cc_ylabels = c(0.52, 0.54, 0.56, 0.58, 0.60, 0.62, 0.64, 0.66)
    ny_ylabels = c(0.30, 0.32, 0.34, 0.36, 0.38, 0.40, 0.42)
  }
  else if (outcomes[o]=='Any Auto Loan or Lease') {
    cc_ylabels = c(0.1, 0.12, 0.14, 0.16, 0.18, 0.2, 0.22)
    ny_ylabels = c(0.11, 0.12, 0.13, 0.14, 0.15, 0.16)
  }
  else if (outcomes[o]=='Total Balance: Col. and Del.') {
    cc_ylabels = c(1800, 2000, 2200, 2400, 2600, 2800, 3000)
    ny_ylabels = c(2000, 2200, 2400, 2600, 2800, 3000)
  }
  else if (outcomes[o]=='Financial Health Index') {
    cc_ylabels = c(-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2)
    ny_ylabels = c(-0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15)
  }
  
  cc_ylimits = c(min(cc_ylabels), max(cc_ylabels))
  ny_ylimits = c(min(ny_ylabels), max(ny_ylabels))
  
  
  g_1 = ggplot(df_cc) + 
    theme_classic() +
    geom_vline(xintercept = 0, color = "black", alpha = .75, linetype = 'dashed') +
    geom_line(aes(x = time, y = evicted,  color = "Evicted", linetype = "Evicted"), size = 2)  +
    geom_line(aes(x = time, y = non_evicted, color = "Not Evicted", linetype = "Not Evicted"), size = 2) +
    scale_color_manual(name = "", 
                       values = c("Evicted" = "#1B9E77", "Not Evicted" = "#D95F02")) +
    scale_linetype_manual(name = "",
                          values = c("Evicted" = "solid", "Not Evicted" = "dashed")) +
    custom_theme + 
    scale_x_continuous(breaks = x_breaks, labels = x_labs) +
    scale_y_continuous(limits = cc_ylimits, breaks = cc_ylabels, labels = cc_ylabels) +
    labs(title='Cook County', x=ifelse(relative_quarter[o], 'Relative Quarter', 'Relative Year') , y=outcomes[o])
  
  g_2 = ggplot(df_ny) + 
    theme_classic() +
    geom_vline(xintercept = 0, color = "black", alpha = .75, linetype = 'dashed') +
    geom_line(aes(x = time, y = evicted,  color = "Evicted", linetype = "Evicted"), size = 2)  +
    geom_line(aes(x = time, y = non_evicted, color = "Not Evicted", linetype = "Not Evicted"), size = 2) +
    scale_color_manual(name = "", 
                       values = c("Evicted" = "#1B9E77", "Not Evicted" = "#D95F02")) +
    scale_linetype_manual(name = "",
                          values = c("Evicted" = "solid", "Not Evicted" = "dashed")) +
    
    custom_theme + 
    scale_x_continuous(breaks = x_breaks, labels = x_labs) +
    scale_y_continuous(limits = ny_ylimits, breaks = ny_ylabels, labels = ny_ylabels) +
    labs(title='New York', x=ifelse(relative_quarter[o], 'Relative Quarter', 'Relative Year'), y='')
  
  g = ggarrange(g_1, g_2, nrow=1, ncol=2, common.legend=TRUE, legend="bottom")
  
  all_graphs[[outcomes[o]]] = g
  
  ggsave(paste0(save_dir, save_names[o], ".png"), plot=g, width=15, height=6)
}





################################################################################
#
## SUPPLEMENTARY FIGURES F.1 and F.2
# 
################################################################################

# Graph-specific settings
boxes = c("earn_CC", "empl_CC", "earn_NY", "empl_NY")
boxes_labels = c("Earnings", "Employment","Earnings", "Employment")
file_names = c('Supplementary data_Figure_F1_Earnings_CC', 'Supplementary data_Figure_F2_Employment_CC', 
               'Supplementary data_Figure_F1_Earnings_NY', 'Supplementary data_Figure_F2_Employment_NY')
plot_digs = c(0,3,0,3)
num_ticks = c(6,5,6,5)

for (o in 1:length(boxes)){
  
  print(boxes[o])
  
  df = boxplots %>% filter(outcome==boxes[o])
  
  if (min(df$beta) < 0){
    
    ub = 0
    lb = min(df$beta) - max(df$se)
  }
  
  if (min(df$beta) > 0){
    ub = max(df$beta) + max(df$se)
    lb = 0
  }
  
  if (boxes[o]=="earn_CC" | boxes[o]=="empl_CC" ) {
    limit_labels = c("Evicted vs ACS", "Evicted vs Non-Evicted", "DiD")
    fill_colors = c("grey25", "grey55", "grey85")
  }
  
  if (boxes[o]=="earn_NY" | boxes[o]=="empl_NY" ) {
    limit_labels = c("Evicted vs Benefits", "Evicted vs Non-Evicted", "DiD")
    fill_colors = c("grey25", "grey55", "grey85")
  }
  
  z = ggplot(df) +
    theme_classic() + 
    geom_hline(aes(yintercept=0), color="black") +
    geom_bar( aes(x=id, y=beta), stat="identity", 
              fill= fill_colors, color="black",
              alpha=0.8) +
    geom_errorbar( aes(x=id, y=beta, ymin=beta-se, ymax=beta+se), 
                   width=0.2, colour="black", alpha=0.9, size=1) +
    custom_theme + 
    theme(strip.text.x = element_text(size = 25)) +
    theme(strip.background = element_blank()) + 
    scale_y_continuous(breaks = waiver(), n.breaks=5) +
    scale_x_discrete(labels = limit_labels) + 
    labs(x='', y=boxes_labels[o])
  
  ggsave(paste0(save_dir, file_names[o], ".png"), plot=z, width=11, height=6)
  
}

